classdef finiteDifferentiator
  %FINITEDIFFERENTIATOR Flexible class for finite difference approximation
  %   Takes a (1-dimensional) grid of values, computes finite difference
  %   coefficients once in an initialization step, actual execution should
  %   then be fast
  
  properties (GetAccess = public, SetAccess = immutable)
    xGrid; % grid of x values
    nApproxPoints; %number of approximation points for sufficiently inner point, assumed to be odd
    d1Matrix; % matrix for first difference formula
    d2Matrix; % matrix for second difference formula
  end
  
  methods (Access = public)
    function obj = finiteDifferentiator(xGrid,nApproxPoints)
      %FINITEDIFFERENTIATOR Construct an instance of this class
      %   Checks inputs, computes the finite difference formulas and stores
      %   everything in the properties for later reference
      %   assumes that length(xGrid)>=nApproxPoints
      if nApproxPoints < 3 || mod(nApproxPoints,2)~=1
        error('Input ''nApproxPoints'' is %d, but odd integer >=3 was expected!',nApproxPoints);
      end
      obj.nApproxPoints = nApproxPoints;
      if size(xGrid,1)==1
        xGrid = xGrid';
      end
      obj.xGrid = xGrid;
      nNeighborsUsed = (nApproxPoints-1)/2;
      % add artificial points to boundaries to have sufficient neighbors in
      % all cases
      xGridLong = [xGrid(1)*ones(nNeighborsUsed,1);xGrid;xGrid(end)*ones(nNeighborsUsed,1)];
      xStartIdx = nNeighborsUsed+1;
      xEndIdx = numel(xGridLong) - nNeighborsUsed;
      d1Matrix = spalloc(numel(xGrid),numel(xGrid),numel(xGrid)*nApproxPoints);
      d2Matrix = spalloc(numel(xGrid),numel(xGrid),numel(xGrid)*nApproxPoints);
      neighborIndexOffset = -nNeighborsUsed:nNeighborsUsed;
      for i=xStartIdx:xEndIdx
        % set up equation system matrix (same for first and second
        % difference formula)
        systemMatrix = NaN(nApproxPoints,nApproxPoints);
        for j=1:nApproxPoints
          termOrder = j-1; %order of term this equation row represents
          if termOrder <= nNeighborsUsed+i-xStartIdx && termOrder <= nNeighborsUsed+xEndIdx-i
            systemMatrix(j,:) = (xGridLong(i+neighborIndexOffset)' - xGridLong(i)).^termOrder;
          else
            systemMatrix(j,:) = zeros(1,nApproxPoints);
            if termOrder>nNeighborsUsed+i-xStartIdx
              systemMatrix(j,termOrder-(nNeighborsUsed+i-xStartIdx)) = 1;
            else
              systemMatrix(j,nApproxPoints+1 + nNeighborsUsed+xEndIdx-i - termOrder) = 1;
            end
          end
        end
        rhsFirstDifferences = [0;1;zeros(nApproxPoints-2,1)];
        rhsSecondDifferences = [0;0;2;zeros(nApproxPoints-3,1)];
        solution1d = systemMatrix\rhsFirstDifferences;
        solution2d = systemMatrix\rhsSecondDifferences;
        coeffIdxLeft = 1-min(0,i-nNeighborsUsed-xStartIdx);
        coeffIdxRight = nApproxPoints-max(0,i+nNeighborsUsed-xEndIdx);
        rowIndices = i+neighborIndexOffset-xStartIdx+1;
        rowIndices = rowIndices(rowIndices > 0 & rowIndices <= numel(xGrid));
        d1Matrix(i-xStartIdx+1,rowIndices) = solution1d(coeffIdxLeft:coeffIdxRight)';
        d2Matrix(i-xStartIdx+1,rowIndices) = solution2d(coeffIdxLeft:coeffIdxRight)';
      end
      obj.d1Matrix = d1Matrix;
      obj.d2Matrix = d2Matrix;
    end
    
    function dy = d1(obj,y)
      %D1 computes first derivative of y
      %   y is assumed to be a column vector of same dimension as xGrid
      dy = obj.d1Matrix*y;
    end
    
    function d2y = d2(obj,y)
      %D2 computes second derivative of y
      %   y is assumed to be a column vector of same dimension as xGrid
      d2y = obj.d2Matrix*y;
    end
  end
end

